1 Objetivo

Comparar munícipios para estabelecer metas de redução de óbito e acidentes em estradas, para o ano de 2017. A técnica de clusterização foi utilizada para separar grupos com variáveis mais relevantes umas com as outras.

2 Variáveis do dataset

2.0.1 Carregando os pacotes

library(abjutils)
library(brazilmaps)
library(corrgram)
library(DataExplorer)
library(dplyr)
library(DT)
library(ggplot2)
library(magrittr)
library(maps)
library(mice)
library(plotly)
library(randomForest)
library(readr)
library(readxl)
library(stringr)
library(viridis)

2.0.2 Leitura dos dados

df_dataset <- read_xlsx("city_dataset.xlsx")
df_municipios <- read_xls("MunicipiosBrasil.xls")

2.0.3 Tabela com as variáveis do dataset

df_table <- df_dataset %>% dplyr::select(-ibge)

sketch = htmltools::withTags(table(
 tableHeader(c('Município', 'Ano','PIB', 'Mat1517', 'Veículos', 'Motos',
               'População', 'Pop1519', 'Pop2024', 'Pop2529', 'Pop60p',
               'jovem', 'pjovem', 'pmotos', 'pmat', 'rodovia'))
))

DT::datatable(container = sketch, df_table, rownames = FALSE,
                        options = list(searching = TRUE, pageLength = 10,
                        lengthMenu = c(20,50,100), dom = 'Bfrtip',
                         
                          initComplete = JS("
function(settings, json) {
$(this.api().table().header()).css({
'background-color': '#000',
'color': '#fff'
});
}")))
  • Cabe mencionar que as variáveis com proporção maiores que 1 foram excluídas da análise.

2.0.4 Leitura de base externa

Dados do IBGE foram utilizados para obtenção das variáveis latitude e longitude para a elaboração de um mapa exploratório simples.

df_dataset %<>% rename(municipio = cidade)
df_municipios %<>% rename(municipio = MUNICIPIO)

#Colocando os dados em formato minúsculo
df_municipios$municipio %<>% str_to_title()
#Removendo acentos
df_dataset$municipio %<>% 
    str_to_lower() %>% #transformando em letras minúsculas para facilitar a remoção dos acentos
    str_replace_all("á", "a") %>% 
    str_replace_all("é", "e") %>% 
    str_replace_all("í", "i") %>% 
    str_replace_all("ó", "o") %>% 
    str_replace_all("ú", "u") %>% 
    str_replace_all("ã", "a") %>% 
    str_replace_all("â", "a") %>% 
    str_replace_all("ê", "e") %>% 
    str_replace_all("ô", "o") %>% 
    str_replace_all("õ", "o") %>% 
    str_replace_all("ç", "c") %>% 
    str_to_title() #voltando com as variáveis iniciando com letra maiúscula

df_map <- df_dataset %>% 
  dplyr::select(municipio, populacao, ano) %>% 
  dplyr::distinct() %>% 
  left_join(df_municipios, by = c("municipio"))

df_map %<>% rename(long = LONGITUDE, lat = LATITUDE, pop = populacao)
df_map %<>% dplyr::filter(ano == 2014)

3 Criação de um pequeno mapa com as coordenadas de Latitude e Longitude

p2 <-  df_map %>%
  dplyr::filter(municipio != "Sao Paulo") %>% 
  arrange(pop) %>%
  mutate( name=factor(municipio, unique(municipio))) %>%
  mutate( mytext=paste("Município: ", municipio, "\n", "População: ", pop, sep="")) %>% 
  ggplot() +
    geom_polygon(aes(x=long, y = lat), fill="grey", alpha=0.3) +
    geom_point(aes(x=long, y=lat, size=pop, color = pop, text=mytext, alpha=pop)) +
    scale_size_continuous(range=c(1,15)) +
    scale_color_viridis(option="inferno", trans="log") +
    scale_alpha_continuous(trans="log") +
    theme_void() +
    coord_map() +
    labs(colour = "População")

ggplotly(p2, tooltip = "text", width = 600, height = 800)

4 Imputação por Random Forest

Na técnica de imputação por Random Forest são combinados resultados de diversos classificadores, sendo o resultado um emsemble de várias árvores de decisão. Cada classificador funciona como um preditor.

PS: Para que o documento não ficasse com mais de 1500 linhas de imputação, o output não foi mostrado.

#set.seed(2709) #configurando uma semente
#df_dataset %<>%  filter(pmat <= 1)
#imp <- mice(df_dataset, m = 5, maxit = 50, method = "rf", seed = 2709)
#
#df_dados_imputados <- mice::complete(imp, 1)
#df_dados_imputados %<>% dplyr::filter(ano == 2016)
#Salvando a base para que não seja rodada novamente pelo alto custo computacional
#df_dados_imputados %>% write.csv2("dados_imputados.csv", row.names = FALSE)
df_dados_imputados <- read.csv2("dados_imputados.csv")
df_dados_imputados %<>%  filter(pmat <= 1)

5 Correlação entre as variáveis

Conforme o esperado, as variáveis de população são perfeitamente correlacionadas entre si.

df_dados_imputados %>% 
  dplyr::select(-municipio, -ano, -ibge) %>% 
  corrgram(order=TRUE, upper.panel=panel.cor, main="Gráfico de correlação")

6 Técnica de Análise de Componentes Principais (ACP)

A técnica de ACP constiste em transformar um conjunto de variáveis em diversas componentes principais. Cada componente é uma combinação linear das variáveis originais, sendo independentes entre si. A ideia é reduzir o máximo de informações, com uma perda ínfima de informação.

Na análise foram selecionadas 5 componentes que representam 98.7% da variabilidade dos dados, o que pode ser visto no gráfico abaixo.

label <- as.factor(df_dados_imputados[[1]])

covtrain <- df_dados_imputados %>% 
            dplyr::select(-municipio, -ano, -ibge) %>% 
            scale() %>% 
            cov()

train_pc <- prcomp(covtrain)
varex <- train_pc$sdev^2/sum(train_pc$sdev^2)
varcum <- cumsum(varex)
result <- data.frame(num=1:length(train_pc$sdev),
                         ex=varex,
                         cum=varcum)

plot(result$num,result$cum,type="b",xlim=c(0,14),
     main="Variância explicada pelas 14 componentes",
     xlab="Número de Componentes",ylab="Variância explicada")
abline(v = 5, lty = 2)

Considerando o município de São Paulo na análise. * No gráfico abaixo são plotadas as duas primeiras componentes principais que representam 88% dos dados.

train_score <-  df_dados_imputados %>%
                dplyr::select(-municipio, -ano, -ibge) %>% 
                as.matrix() %>% 
                scale() %*% train_pc$rotation[,1:2]
            
train <- cbind(label,as.data.frame(train_score))

colors <- rainbow(length(unique(train$label)))
names(colors) <- unique(train$label)
plot(train$PC1,train$PC2,type="n",main="Duas primeiras componentes principais",
     xlab = "Segunda componente principal", ylab = "Primeira componente principal")
text(train$PC1,train$PC2,label=train$label, col=colors[train$label])

rm(label)
rm(covtrain)
rm(train_pc)
rm(varex)
rm(varcum)
rm(result)

7 Outlier: São Paulo

Com 6 componentes explicamos aproximadamente 99.5% da variabilidade dos dados.

df_save_data <- df_dados_imputados
df_dados_imputados %<>% dplyr::filter(municipio != "Sao Paulo")
label <- as.factor(df_dados_imputados[[1]])

covtrain <- df_dados_imputados %>% 
            dplyr::select(-municipio, -ano, -ibge) %>% 
            scale() %>% 
            cov()

train_pc <- prcomp(covtrain)
varex <- train_pc$sdev^2/sum(train_pc$sdev^2)
varcum <- cumsum(varex)
result <- data.frame(num=1:length(train_pc$sdev),
                         ex=varex,
                         cum=varcum)

plot(result$num,result$cum,type="b",xlim=c(0,14),
     main="Variância explicada pelas 14 componentes",
     xlab="Número de Componentes",ylab="Variância explicada")
abline(v = 6, lty = 2)

No gráfico de componentes abaixo, sem a presença do município de SP, verificamos um melhor espalhamento dos demais municípíos.

Municípios com maior PIB e maior população como Campinas, Guarulhos, Ribeirão Preto e São Bernardo do Campo tendem a se concentrar mais na primeira componente principal.

Enquanto municípios menores em população e PIB como Pracinha, Reginópolis e Lavínia se concentram na segunda componente principal

train_score <-  df_dados_imputados %>%
                dplyr::select(-municipio, -ano, -ibge) %>% 
                as.matrix() %>% 
                scale() %*% train_pc$rotation[,1:2]
            
train <- cbind(label,as.data.frame(train_score))

colors <- rainbow(length(unique(train$label)))
names(colors) <- unique(train$label)
plot(train$PC1,train$PC2,type="n",main="Duas primeiras componentes principais",
     xlab = "Segunda componente principal", ylab = "Primeira componente principal")
text(train$PC1,train$PC2,label=train$label, col=colors[train$label])

8 Modelo de regressão - Variáveis relevantes

set.seed(2809)
df_lm <- df_dados_imputados %>% 
  na.omit() %>% 
  dplyr::select(-municipio, -ano, -ibge) %>% 
  mutate(pib = log(pib), mat1517 = log(mat1517), veiculos = log(veiculos), 
         motos = log(motos), populacao = log(populacao), pop1519 = log(pop1519),
         pop2024 = log(pop2024), pop2529 = log(pop2529), jovem = log(jovem))

regressao <- lm(log(pib) ~ ., data = df_lm %>% 
                       dplyr::select(-populacao,-pop1519, -jovem, -pop2529,-pop2024, -motos, -pmotos))
summary(regressao)
## 
## Call:
## lm(formula = log(pib) ~ ., data = df_lm %>% dplyr::select(-populacao, 
##     -pop1519, -jovem, -pop2529, -pop2024, -motos, -pmotos))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.215678 -0.028185 -0.001953  0.028722  0.217919 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.117e+00  8.911e-02  23.751  < 2e-16 ***
## mat1517      4.673e-02  1.699e-02   2.750 0.006201 ** 
## veiculos     4.467e-02  3.677e-03  12.150  < 2e-16 ***
## pop60p      -4.617e-03  1.206e-03  -3.827 0.000148 ***
## pjovem      -2.857e-03  1.415e-03  -2.019 0.044069 *  
## pmat        -8.617e-02  1.615e-02  -5.337 1.51e-07 ***
## rodovia      2.552e-04  7.602e-05   3.357 0.000856 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04675 on 443 degrees of freedom
## Multiple R-squared:  0.8007, Adjusted R-squared:  0.798 
## F-statistic: 296.7 on 6 and 443 DF,  p-value: < 2.2e-16

8.0.1 Multicolinearidade

  • A variável jovem que representa a população de 15 a 19, a de 20 e 24 e a 25 e 29 anos, foi retirada do modelo, por apresentar multicolinearidade. O que é confirmada com a análise do VIF logo abaixo.

  • O VIF (fator de variância generalizada) avalia se há multicolinearidade entre as variáveis, isto é, se as variáveis são altamente correlacionadas, o que aumenta a variabilidade (desvio padrão) do modelo.

  • Se o VIF < 5, a variável permanece no modelo.

  • Logo as variáveis selecionadas explicam bem o modelo.eu

car::vif(regressao)
##  mat1517 veiculos   pop60p   pjovem     pmat  rodovia 
## 1.499978 4.844512 2.441433 2.265945 3.992041 1.416054

9 Clusterização - Kmeans clustering

A ideia da técnica é encontrar agrupamentos nos dados. O algoritmo trabalha iterativamente para ‘encaixar’ cada observação em um cluster baseado nas características que são providenciadas no dataset.

A distância entre os agrupamentos é a distância entre os centróides que são os valores médios das observações. Toda vez que um agrupamento é formado, um novo método é calculado.

Cada ponto é assinalado ao seu centróide mais próximo, baseado no quadrado da distância Euclideana.

Para que a análise não seja distorcida, as variáveis foram padronizadas, a fim de tirar a medida em grandezas diferentes.

Logo após a padronização, foi rodado o algoritmo K-Means.

9.0.1 Número de clusters

Abaixo é feita um código para obter a soma dos quadrados dentro dos clusteres. Uma forma de se escolher o número ideal de clusters é a partir que a a soma dos começa a estabilizar, o que acontece com 7 ou 8 clusters.

set.seed(2709) #configurando a semente
#Voltando com os dados originais, incluindo São Paulo
df_dados_imputados <- df_save_data

# Padronizando os dados
df_scale <- df_dados_imputados %>%
                dplyr::select(-municipio, -ano, -ibge) %>% 
                scale() #Distância Euclideana

wss <- (nrow(df_scale)-1)*sum(apply(df_scale,2,var))

for (i in 2:15) wss[i] <- sum(kmeans(df_scale,
   centers=i)$withinss)

plot(1:15, wss, type="b", xlab="Número de clusters",
  ylab="Soma dos quadrados intra cluster", main = "Escolha do número clusters")
abline(v = 6, col = "red", lty = c(2), lwd = c(3))

cluster_km <- kmeans(df_scale, centers = 6)$cluster

9.0.2 Número de municípios/Proporção de municípios por cluster

df_dados_imputados %>%
  dplyr::select(municipio, -ano, -ibge) %>%
  cbind(cluster_km) %>%
  rename(Cluster = cluster_km) %>% 
  arrange(Cluster) %>%
  mutate(N = 1) %>% 
  group_by(Cluster) %>% 
  summarize(Total_Municipios = sum(N)) %>%
  mutate(Proporcao_Municipios = round(100*Total_Municipios/sum(Total_Municipios),2))
## # A tibble: 6 x 3
##   Cluster Total_Municipios Proporcao_Municipios
##     <int>            <dbl>                <dbl>
## 1       1               1.                0.220
## 2       2             127.               28.2  
## 3       3             129.               28.6  
## 4       4               8.                1.77 
## 5       5              42.                9.31 
## 6       6             144.               31.9

9.0.3 Tabela com a média de cada variável para cada cluster

df_dados_imputados %>%
  dplyr::select(-municipio, -ano, -ibge) %>%
  cbind(cluster_km) %>% 
  rename(Cluster = cluster_km) %>% 
  group_by(Cluster) %>% 
  summarize_all(mean)
## # A tibble: 6 x 15
##   Cluster    pib mat1517 veiculos  motos populacao pop1519 pop2024 pop2529
##     <int>  <dbl>   <dbl>    <dbl>  <dbl>     <dbl>   <dbl>   <dbl>   <dbl>
## 1       1 5.02e8    79.2 7010508. 9.65e5 11638802. 845272. 874967. 975688.
## 2       2 1.73e6    68.4   31979. 6.84e3    66953.   5598.   5767.   5972.
## 3       3 2.20e6    77.6   44070. 1.18e4    64965.   4763.   5113.   5376.
## 4       4 3.13e5    58.2    2728. 7.41e2     8610.    607.   1165.   1173.
## 5       5 9.84e6    79.1  163508. 3.45e4   291973.  21676.  23209.  25018.
## 6       6 2.41e5    63.6    5580. 1.13e3    11886.    894.    961.    977.
## # ... with 6 more variables: pop60p <dbl>, jovem <dbl>, pjovem <dbl>,
## #   pmotos <dbl>, pmat <dbl>, rodovia <dbl>

9.0.4 Análise dos clusters

  • Cluster 1: São Paulo.

O cluster 1 tem apenas São Paulo como município. Por apresentar as maiores médias do PIB e do número de veículos, maior número de rodovias em KM, a menor proporção de matriculados no Ensino Médio na faixa de 15 a 19 anos e a menor proporção da população de 60 anos ou mais e a maior proporção da população na faixa de idade entre 25 a 29 anos, aparece sozinho no agrupamento.

  • Cluster 2: Boituva, Cabreúva, Cubatão, Hortolândia, Taboão da Serra.

Tem um percentual baixo da população abaixo de 60 anos e um percentual alto da população jovem.

  • Cluster 3: Catanduva, Mirassol, Mogi Mirim, Ourinhos, Ubatuba.

Tem um PIB médio pequeno e um percentual de população entre 15 e 17 anos alto.

  • Cluster 4: Balbinos, Lavínia, Reginópolis

Tem um número baixo tanto de veículos como motos, uma população pequena comparada aos demais municípios e a que tem a maior proporção jovens na população.

  • Cluster 5: Araraquara, Guarulhos, Barretos, Itapeva, Franca, Jundiaí.

O segundo cluster que tem o maior número de jovens, perdendo apenas para o cluster 1 (município de São Paulo); a menor proporção de jovens matriculados no ensino médio, só perdendo para São Paulo;

  • Cluster 6: Cardoso, Dourado, Divinolândia, Nova Granada, Vista Alegre do Alto.

São municípios pequenos em população, PIB médio baixo, grande parte da população na faixa dos 60 anos de idade ou mais.